import sys
sys.path.append('C:/Users/Hu/Dropbox/Research/PythonWork/Cancer/src/STAT/')

import numpy as np
import csv
from datetime import datetime
from time import strftime
import matplotlib.pyplot as plt
import math
import random
from scipy import stats
import ols
from scipy.stats import f


def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def build_data_list(inputCSV, datatype="int"):
    sKey = []
    fn = inputCSV
    f = open(inputCSV)
    #ra = csv.DictReader(file(fn), dialect="excel")
    ra = csv.DictReader(f, dialect="excel")

    if datatype=="int":
        for record in ra:
            #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
            for item in ra.fieldnames:
                temp = int(float(record[item]))
                sKey.append(temp)
    elif datatype=="float":
        for record in ra:
            for item in ra.fieldnames:
                temp = float(record[item])
                sKey.append(temp)
    else:
        print 'incorrect value in second parameter (datatype)'
        return


    sKey = np.array(sKey)
    sKey.shape=(-1,len(ra.fieldnames))
    return sKey

def geneBpChoices(tmin, tmax, num):
    tmin = float(tmin)
    tmax = float(tmax)
    num = int(num)
    step = (tmax - tmin)/num
    tbpchoices = []
    for i in range(num):
        tbpchoices.append(tmin+i*step)
    return tbpchoices

def pickBreakpoint(response, predictor, bpvarible, bpChoices, predictorName):
    results = np.zeros((len(bpChoices)-1, 2))
    print bpChoices
    bpvarible = np.array(bpvarible)
    predictor = np.array(predictor)
    predictor_t = np.transpose(predictor)
    predictorName = np.append(predictorName, 'bpvarible')
    #print predictor_t, predictor_t.shape
    
    for i in range(len(bpChoices)-1):
        x2star = (bpvarible - bpChoices[i]) * np.greater(bpvarible, bpChoices[i])
        tempPredictor = np.append(predictor_t, x2star)
        tempPredictor.shape = (len(predictorName),-1)
        tempPredictor = np.transpose(tempPredictor)
        #print tempPredictor
        tempmodel = ols.ols(response, tempPredictor,'y', predictorName)
        results[i,0] = i
        results[i,1] = tempmodel.R2
    print results

    optBP = int(results[np.argmax(results, axis = 0)[1],0])
    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', results[optBP, 1]

    x2star = (bpvarible - bpChoices[optBP]) * np.greater(bpvarible, bpChoices[optBP])
    tempPredictor = np.append(predictor_t, x2star)
    tempPredictor.shape = (len(predictorName),-1)
    optPredictor = np.transpose(tempPredictor)
    optmodel = ols.ols(response, optPredictor,'y', predictorName)
    #optmodel = ols.ols(response, optPredictor,'y',predictorName)
    
    y_hat = np.dot(optmodel.x, optmodel.b)
    #fileLoc = filepath + 'pieceR_model2_y_hat.csv'
    #np.savetxt(fileLoc, y_hat, delimiter=',', fmt = '%s')
    
    print optmodel.summary()
    print 'MSE =', optmodel.sse
    #print 'before bp, b0 =',optmodel.b[0], ', b1 =', optmodel.b[1], ', bd =', optmodel.b[2]
    #print 'after bp, b0 =', optmodel.b[0] - optmodel.b[3] * bpChoices[optBP], ', b1 =', optmodel.b[1], ', bd =', optmodel.b[2] + optmodel.b[3]

    #calpredictedvalue(zip(x1, predictor), bpChoices[optBP], zip(optmodel.b, [optmodel.b[0] - optmodel.b[3] * bpChoices[optBP], optmodel.b[1], optmodel.b[2] + optmodel.b[3]]), 'exp_inoutflow_model2A.csv')
    return results, y_hat


    results = np.zeros((len(bpChoices)-1, 2))
    print bpChoices
    
    for i in range(len(bpChoices)-1):
        print i
        x2star = (predictor - bpChoices[i]) * np.greater(predictor, bpChoices[i])
        x1star = x1 * np.greater(predictor, bpChoices[i]) 
        tempPredictor = np.array(zip(x1, x1star, predictor, x2star))
        #fileLoc = filePath + 'temp.csv'
        #np.savetxt(fileLoc, tempPredictor, delimiter=',', fmt = '%s')
        #print tempPredictor
        tempmodel = ols.ols(response, tempPredictor,'y' ,predictorName)
        results[i,0] = i
        #results[i,1] = tempmodel.sse
        results[i,1] = tempmodel.R2

    optBP = int(results[np.argmax(results, axis = 0)[1],0])
    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', results[optBP, 1]

    #x2star = (predictor - bpChoices[optBP]) * np.greater(predictor, bpChoices[optBP])
    #optPredictor = np.array(zip(predictor, x2star))
    #optmodel = ols.ols(response, optPredictor,'y',['x1', 'x2'])
    x1star = x1 * np.greater(predictor, bpChoices[optBP])
    x2star = (predictor - bpChoices[optBP]) * np.greater(predictor, bpChoices[optBP])
    optPredictor = np.array(zip(x1, x1star, predictor, x2star))
    optmodel = ols.ols(response, optPredictor,'y',predictorName)
    
    #return bpChoices[optBP], results, optmodel, optmodel.b[0]+optmodel.b[1]*predictor+optmodel.b[2]*x2star
    print results, optmodel.b
    print optmodel.summary()
    return results

def separateData(data1, data2, dataC, threshold):
    if (len(data1) <> len(data2)) or (len(data1) <> len(dataC)):
        print 'Error in input data length! --- seperateData.error.1'

    data1Left = []
    data1Right = []
    data2Left = []
    data2Right = []
    dataleftIdx = []
    datarightIdx = []
    
    for i in range(len(data1)):
        if dataC[i] < threshold:
            data1Left.append(data1[i])
            data2Left.append(data2[i])
            dataleftIdx.append(i)
        else:
            data1Right.append(data1[i])
            data2Right.append(data2[i])
            datarightIdx.append(i)

    return np.array(data1Left), np.array(data1Right), np.array(data2Left), np.array(data2Right), dataleftIdx, datarightIdx

def calculateR2(tleftmodel, trightmodel, totalmean):
    #self.sse = dot(self.e,self.e)/self.df_e
    #tsse = np.dot(tmodel.e, tmodel.e)/tmodel.df_e
    #tsse = np.dot(tmodel.e, tmodel.e)
    tsse = tleftmodel.sse * tleftmodel.df_e + trightmodel.sse * trightmodel.df_e
    
    #y_hat = np.dot(tmodel.x,tmodel.b)
    tsst = np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean)
    tsst = tsst + np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    
    #print tleftmodel.sse * tleftmodel.df_e, trightmodel.sse * trightmodel.df_e, np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean), np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    #print 1 - (tleftmodel.sse * tleftmodel.df_e)/np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean)
    #print 1 - (trightmodel.sse * trightmodel.df_e)/np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    
    #return 1-tsse/tsst
    #print ((tleftmodel.nobs+trightmodel.nobs-1.0), (tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef)), ((tleftmodel.nobs+trightmodel.nobs-1.0)/(tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef))
    return 1 - (tsse/tsst)*((tleftmodel.nobs+trightmodel.nobs-1.0)/(tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef))

def simpleRegression(response, predictor, predictorName):
    model = ols.ols(response, predictor,'y',predictorName)
    #model = ols.ols(response, predictor)
    print model.summary()
    print model.b
    print predictor
    #print 'b0 =',model.b[0], ', b1 =', model.b[1], ', bd =', model.b[2], ', MSE =', model.sse
    return np.dot(model.x, model.b), model

def separateRegression(response, predictor, sepData, bpChoices, predictorName, equaltestid, equaltestid2):
    results = np.zeros((len(bpChoices)-1, 4))
    print bpChoices
    chosencriterion = 'r2'  # max(r2) is the same criterion as min(sse)
    for bpid in range(len(bpChoices)-1):
        #print bpid
        responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(response, predictor, sepData, bpChoices[bpid])
        #print predictorLeft
        leftmodel = ols.ols(responseLeft, predictorLeft,'y',predictorName)
        rightmodel = ols.ols(responseRight, predictorRight,'y',predictorName)
        results[bpid,0] = bpid
        if chosencriterion == 'r2':           
            results[bpid,1] = leftmodel.R2adj
            results[bpid,2] = rightmodel.R2adj
            #results[bpid,3] = 1 - (leftmodel.e.var() + rightmodel.e.var())/(leftmodel.y.var() + rightmodel.y.var())
            results[bpid,3] = calculateR2(leftmodel, rightmodel, np.mean(response))
            #results[bpid,3] = (leftmodel.R2 + rightmodel.R2)/2
        elif chosencriterion == 'sse':
            results[bpid,1] = leftmodel.sse * leftmodel.df_e
            results[bpid,2] = rightmodel.sse * rightmodel.df_e
            results[bpid,3] = results[bpid,1] + results[bpid,2]
            #yhatL = np.dot(leftmodel.x, leftmodel.b)
            #print results[bpid,1], np.sum((responseLeft - yhatL)**2)

    print results
    if chosencriterion == 'r2':
        optBP = int(results[np.argmax(results, axis = 0)[-1],0])
    elif chosencriterion == 'sse':
        optBP = int(results[np.argmin(results, axis = 0)[-1],0])

    responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(response, predictor, sepData, bpChoices[optBP])
    leftmodel = ols.ols(responseLeft, predictorLeft,'y',predictorName)
    rightmodel = ols.ols(responseRight, predictorRight,'y',predictorName)

    #equaltestid = int(equaltestid)
    if equaltestid[0] > -1:
        temppredictorLeft = predictorLeft
        temppredictorLeft[:,equaltestid[1]] = temppredictorLeft[:,equaltestid[0]] + temppredictorLeft[:,equaltestid[1]]
        temppredictorLeft = np.delete(temppredictorLeft, equaltestid[0], 1)
        temppredictorRight = predictorRight
        temppredictorRight[:,equaltestid[1]] = temppredictorRight[:,equaltestid[0]] + temppredictorRight[:,equaltestid[1]]
        temppredictorRight = np.delete(temppredictorRight, equaltestid[0], 1)
        temppredictorName = np.delete(predictorName, equaltestid[0], None)
        templeftmodel = ols.ols(responseLeft, temppredictorLeft,'y',temppredictorName)
        temprightmodel = ols.ols(responseRight, temppredictorRight,'y',temppredictorName)
        fleft = (leftmodel.R2 - templeftmodel.R2) * (leftmodel.nobs - len(predictorName) - 2) / (1 - leftmodel.R2)
        fright = (rightmodel.R2 - temprightmodel.R2) * (rightmodel.nobs - len(predictorName) - 2) / (1 - rightmodel.R2)
        pleft = 1 - f.cdf(fleft, 1, leftmodel.nobs - len(predictorName) - 2)
        pright = 1 - f.cdf(fright, 1, rightmodel.nobs - len(predictorName) - 2)

        if equaltestid2[0] > -1:
            temppredictorLeft = predictorLeft
            temppredictorLeft[:,equaltestid2[1]] = temppredictorLeft[:,equaltestid2[0]] + temppredictorLeft[:,equaltestid2[1]]
            temppredictorLeft = np.delete(temppredictorLeft, equaltestid2[0], 1)
            temppredictorRight = predictorRight
            temppredictorRight[:,equaltestid2[1]] = temppredictorRight[:,equaltestid2[0]] + temppredictorRight[:,equaltestid2[1]]
            temppredictorRight = np.delete(temppredictorRight, equaltestid2[0], 1)
            temppredictorName = np.delete(predictorName, equaltestid2[0], None)
            templeftmodel = ols.ols(responseLeft, temppredictorLeft,'y',temppredictorName)
            temprightmodel = ols.ols(responseRight, temppredictorRight,'y',temppredictorName)
            fleft = (leftmodel.R2 - templeftmodel.R2) * (leftmodel.nobs - len(predictorName) - 2) / (1 - leftmodel.R2)
            fright = (rightmodel.R2 - temprightmodel.R2) * (rightmodel.nobs - len(predictorName) - 2) / (1 - rightmodel.R2)
            pleft2 = 1 - f.cdf(fleft, 1, leftmodel.nobs - len(predictorName) - 2)
            pright2 = 1 - f.cdf(fright, 1, rightmodel.nobs - len(predictorName) - 2)

    yhatL = np.dot(leftmodel.x, leftmodel.b)
    yhatR = np.dot(rightmodel.x, rightmodel.b)
    yhat = np.zeros(len(response))

    for i in range(len(yhatL)):
        yhat[dataleftIdx[i]] = yhatL[i]

    for i in range(len(yhatR)):
        yhat[datarightIdx[i]] = yhatR[i]

    yhat = np.exp(yhat)
    fileLoc = filepath + 'separateR_model2_y_hat.csv'
    #np.savetxt(fileLoc, yhat, delimiter=',', fmt = '%s')
    #print 'Optimal Index:', optBP
    #print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', calculateR2(leftmodel, rightmodel, np.mean(response))
   
    print '----------------------------- left model -----------------------------'
    print leftmodel.summary()
    print '----------------------------- right model -----------------------------'
    print rightmodel.summary()

    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with adjR2 = ', results[optBP, -1]
   
    outputstring = 'before bp'
    for i in range(len(predictorName)+1):
        outputstring += ', b' + str(i) + ' = ' + "%.2f" %(leftmodel.b[i]) + '(' + "%.3f" %(leftmodel.se[i])+ ')'
    outputstring += ', with adjR2 = ' + "%.4f" %(leftmodel.R2adj)
    if equaltestid[0] > -1:
        outputstring += ', b1 <> b2 with pvalue = ' + "%.3f" %(pleft)
        if equaltestid2[0] > -1:
            outputstring += ', b3 <> b4 with pvalue = ' + "%.3f" %(pleft2)
    print outputstring

    outputstring = 'after bp'
    for i in range(len(predictorName)+1):
        outputstring += ', b' + str(i) + ' = ' + "%.2f" %(rightmodel.b[i]) + '(' + "%.3f" %(rightmodel.se[i])+ ')'
    outputstring += ', with adjR2 = ' + "%.4f" %(rightmodel.R2adj)
    if equaltestid[0] > -1:
        outputstring += ', b1 <> b2 with pvalue = ' + "%.3f" %(pright)
        if equaltestid2[0] > -1:
            outputstring += ', b3 <> b4 with pvalue = ' + "%.3f" %(pright2)
    print outputstring

    syhat, smodel = simpleRegression(response, predictor, predictorName)

    fsmr = (results[optBP, -1] - smodel.R2) * (smodel.nobs - len(predictorName) - 2) / (1 - results[optBP, -1])
    psmr = 1 - f.cdf(fsmr, 1, smodel.nobs - len(predictorName) - 2)

    print '--------------------Simple Regression Result-----------------------'
    outputstring = ''
    for i in range(len(predictorName)+1):
        outputstring += 'b' + str(i) + ' = ' + "%.2f" %(smodel.b[i]) + '(' + "%.3f" %(smodel.se[i])+ '), '
    outputstring += 'with adjR2 :' + '%.4f' %(smodel.R2adj) + ', separate regression with pvalue = ' + "%.3f" %(psmr)
    print outputstring

    #calpredictedvalue(predictor, bpChoices[optBP], zip(leftmodel.b, rightmodel.b), 'exp_inoutflow_model2B.csv')
    #calconfidenceinterval(predictorLeft, predictorRight, [leftmodel.sse, rightmodel.sse], response, predictor, bpChoices[optBP], zip(leftmodel.b, rightmodel.b), 'ci_model2B.csv')
    return results, yhat

def calconfidenceinterval(leftx, rightx, mse, y, x, bp, b,filename):
    leftx = np.array(zip(np.ones(len(leftx)), leftx[:,0], leftx[:,1]))
    rightx = np.array(zip(np.ones(len(rightx)), rightx[:,0], rightx[:,1]))

    templeftx = np.linalg.inv(np.dot(leftx.T, leftx))
    temprightx = np.linalg.inv(np.dot(rightx.T, rightx))

    b = np.array(b)
    t = [stats.t.ppf(0.975, len(leftx) - len(leftx[0,:])), stats.t.ppf(0.975, len(rightx) - len(rightx[0,:]))]

    yhat = []
    i = 0
    #print len(x)
    for item in x:
        if item[1] < bp:
            bpflag = 0
            tempx = templeftx
        else:
            bpflag = 1
            tempx = temprightx
        temp = b[0, bpflag] + b[1, bpflag] * item[0] + b[2, bpflag] * item[1]
        xh = np.array([1, item[0], item[1]])
        #print np.dot(np.dot(xh, tempx), xh.T), mse[bpflag]*(np.dot(np.dot(xh, tempx), xh.T)
        se = math.pow(mse[bpflag]*(np.dot(np.dot(xh, tempx), xh.T)), 0.5) * t[bpflag]
        pise = math.pow(mse[bpflag]*(1 + np.dot(np.dot(xh, tempx), xh.T)), 0.5) * t[bpflag]
        cil = 0
        cih = 0
        pil = 0
        pih = 0
        if y[i] < (temp-se):
            cil = 1
        if y[i] > (temp+se):
            cih = 1
        if y[i] < (temp-pise):
            pil = 1
        if y[i] > (temp+pise):
            pih = 1
        yhat.append([temp, temp-se, temp+se, temp-pise, temp+pise,cil,cih,pil,pih])
        i += 1

    yhat = np.exp(yhat)
    fileLoc = filepath + filename   #'exp_inoutflow_model2A.csv'
    headerstr = 'yhat,cileft,ciright,pileft,piright,cil,cih,pil,pih'
    #print fileLoc
    np.savetxt(fileLoc, yhat, delimiter=',', header = headerstr, fmt = '%f')
    
    #plt.plot(range(len(y)), y, 'g', label='Observed')
    #plt.plot(range(len(y)), yhat[:,0], 'b', label='Predicted')
    #plt.plot(range(len(y)), yhat[:,1], 'b', label='CI Left')
    #plt.plot(range(len(y)), yhat[:,2], 'b', label='CI Right')
    #plt.plot(range(len(y)), yhat[:,3], 'r', label='PI Left')
    #plt.plot(range(len(y)), yhat[:,4], 'r', label='PI Right')
    #plt.show()
    

def calpredictedvalue(x, bp, b, filename):
    # x: [f1f2, dist]
    # b = [[b00, b10], [b01, b11], [b02, b12]]
    print b
    yhat = []
    b = np.array(b)
    bpflag = 1  # 0: before, 1: after
    for item in x:
        if item[1] < bp:
            bpflag = 0
        else:
            bpflag = 1
        temp = b[0, bpflag] + b[1, bpflag] * item[0] + b[2, bpflag] * item[1]
        #temp = np.dot(item, b[:,bpflag])
        yhat.append(temp)

    yhat = np.array(yhat)
    yhat = np.exp(yhat)
    #print yhat
    #fileLoc = filepath + 'model2b_yhat.csv'
    #np.savetxt(fileLoc, yhat, delimiter=',', fmt = '%f')
    expinoutflow = np.zeros((len(inoutflow), 2)) # [inflow, outflow, obs-exp inflow, obs-exp outflow]
    #print len(yhat), len(data)
    #print zip(data, yhat)
    i = 0
    for item in data:
        expinoutflow[int(item[1]), 1] += yhat[i]
        expinoutflow[int(item[2]), 0] += yhat[i]
        i += 1

    fileLoc = filepath + filename   #'exp_inoutflow_model2A.csv'
    headerstr = 'inflow,outflow'
    #np.savetxt(fileLoc, expinoutflow, delimiter=',', header = headerstr, fmt = '%f')
    
def dataprepare():
    filepath = 'C:/_DATA/migration_census_2000/REDCAPLUS_100Region/'

    #dataCSV = filepath + 'migration_100regs_dist.csv'
    dataCSV = filepath + 'census_migration_100regs_dis.csv'
    data = build_data_list(dataCSV)

    print data

    file = filepath + 'census_migration_100regs_pop.csv'
    #file = filepath + '100reg_pop.csv'
    pop = build_data_list(file)

    file = filepath + '100reg_popdensity.csv'
    tdata = build_data_list(file, "float")
    popden = tdata[:,-1]

    '''
    grossflowCSV = filepath + 'migration_100regs_inoutflow.csv'
    #grossflowCSV = filepath + 'census_migration_100regs_inoutflow.csv'
    grossflow = build_data_list(grossflowCSV)

    inflow = np.zeros(len(grossflow))
    outflow = np.zeros(len(grossflow))

    for item in data:
        outflow[item[1]] += item[-1]
        inflow[item[2]] += item[-1]

    print zip(inflow, outflow)
    '''
        
    output = []
    
    for item in data:
        #dist = int(cal_dist(centroid[item[0], 0], centroid[item[0], 1], centroid[item[1], 0], centroid[item[1], 1]))
        #output.append([item[0], grossflow[item[1], 1], grossflow[item[2], 0], pop[item[1], 1], pop[item[2], 0], item[3]])
        output.append([item[0]/1000, pop[item[1], 1], pop[item[2], 0], popden[item[1]], popden[item[2]], item[3]])

    #print output
    headerstr = 'dist,p1,p2,pd1,pd2,flowVol'
    print headerstr
    np.savetxt(filepath + '100reg_regdata_popden.csv', output, delimiter=',', header = headerstr, fmt = '%s')
    return output
    
#--------------------------------------------------------------------------
#MAIN

if __name__ == "__main__":
    print '===================================================='
    print "begin at " + getCurTime()

    filepath = 'C:/_DATA/migration_census_2000/'
    #file = filepath + '100reg_regressiondata.csv' #[dist,f1,f2,p1,p2,vol]: 1:out, 2:in
    #file = filepath + 'regressiondataeast.csv'
    #file = filepath + 'regressiondatawest.csv'
    file = filepath + 'regressiondata.csv' #[dist,graphdist,f1,f2,p1,p2,vol]: 1:out, 2:in
    data = build_data_list(file)   

    #data = dataprepare()    #[dist,p1,p2,pd1,pd2,flowVol]: 1:out, 2:in
    data = np.array(data)

    lndata = np.log(data)

    dist = data[:,0]
    lndist = lndata[:,0]

    '''
    #file = filepath + 'formatdataeast.csv'
    #file = filepath + 'formatdatawest.csv'
    file = filepath + 'census_migration_100regs.csv'
    formatdata = build_data_list(file)
    
    file = filepath + 'census_migration_100regs_pop.csv'
    #file = filepath + 'census_migration_100regs_inoutflow.csv' #[inflow, outflow]
    massdata = build_data_list(file)
    lnmassdata = np.log(massdata)

    
    output = []
    i = 0
    for item in formatdata:
        oid = item[0]
        did = item[1]
        olnmass = lnmassdata[oid,1]
        dlnmass = lnmassdata[did,0]
        #tlnflow = np.log(item[2])
        output.append([olnmass, dlnmass, olnmass+dlnmass, lndist[i]])
        i += 1
    output = np.array(output)
    '''
    inputPredictor = []
    #i = 0
    for i in range(len(data)):
        inputPredictor.append([lndata[i,4], lndata[i,5], data[i,0]])
    inputPredictor = np.array(inputPredictor)

    #inputPredictor = lndata[:,[3,4,0]]

    bpChoices = range(90,4200,10)
    #bpChoices = np.log(range(90,4200,10))
    #bpChoices = np.log(range(90,2100,10))  #east
    #bpChoices = range(90,2100,10)  #east
    #bpChoices = np.log(range(140,3100,10))  #west
    #bpChoices = range(140,3100,10)  #west

    inputPredictorName = ['lnp1', 'lnp2', 'dist']
    #inputPredictorName = ['lnf1', 'lnf2', 'dist']
    #inputPredictorName = ['lnp1', 'lnp2', 'lndist']
    #inputPredictorName = ['lnf1', 'lnf2', 'lndist']
    #inputPredictorName = ['lnp1p2', 'lndist']
    #yhat, smodel = simpleRegression(lndata[:,-1], inputPredictor, inputPredictorName)
    #inputPredictorName = ['lnf1', 'lnf2', 'lnf12', 'lnf22', 'lndist']
    result, yhat = separateRegression(lndata[:,-1], inputPredictor, dist, bpChoices, inputPredictorName, [0,1], [-1,3])

    #inputPredictorName = ['lnp1', 'lnp2', 'lndist']
    #results, y_hat = pickBreakpoint(lndata[:,-1], inputPredictor, lndist, bpChoices, inputPredictorName)

    #fileLoc = filepath + 'piecewiseR_model2_results.csv'
    #fileLoc = filepath + 'separateR_model2_results.csv'
    headerstr = 'exp'
    fileLoc = filepath + 'temp.csv'
    np.savetxt(fileLoc, yhat, delimiter=',', header = headerstr, fmt = '%s')
    